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Abstract 

A few space agencies are planning missions to the vicinity of the Sun- 
Earth L 2 point, some involving a distributed space system of telescope 
spacecraft, configured in a plane about a hub. An improved understanding 
is developed of the relative motion of such objects in formation flight. The 
telescope equations of motion are written relative to the hub, in terms of 
the hub’s distance from L 2 , and an analytical solution is developed, useful 
for performing orbit control analysis. A halo telescope orbit is investigated, 
with initial conditions selected to avoid resonance excitation. An example 
case of the resulting solution is presented. 

INTRODUCTION 

For centuries, astronomers, mathematicians, and engineers have been interested in the mo- 
tion of objects in the vicinity of the so-called Lagrange libration points. These points, 
equilibrium solutions to the circular restricted three-body problem, provide a starting point 
for analysis using a more sophisticated force model. In particular, Farquhar’s ground break- 
ing work [1] made lucid the results from the previous 200 years, rendering accessible the 
topic of restricted three-body motion. 

Eventually, serious consideration was given to the placement of spacecraft near such 
points. In 1978, the ISEE-3 spacecraft was orbited near the Sun-Earth L\ libration point [2]. 
This mission required an advanced description of the vehicle’s motion, in order to allow 
the development of a control scheme to maintain the spacecraft on its desired trajectory. 
The analysis, as described by Richardson [3], used the classical Lindstedt-Poincare method 
of successive approximations in conjunction with a new technique for obtaining periodic 
solutions to nonlinear systems of differential equations. Richardson’s work provided an 
analytical approximation for some particular periodic solutions of the circular restricted 
problem. The approach introduced a valuable procedure for the removal of secular terms 
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in the context of a halo orbit solution, in which the spacecraft was to maintain a periodic 
trajectory about the libration point. 

In one of Howell’s many publications [4], she presented an overview of various contribu- 
tions to finding periodic solutions to the three-body problem. Her work employed dynamical 
systems theory to support trajectory design. This theory is very beneficial in relating initial 
conditions to a suitable phase-space and to better understanding the geometry of the phase 
space in the three-body problem via stable and unstable manifolds. In [5], there is a sig- 
nificant review of the work to find and describe orbits near collinear libration points. That 
presentation offers a fairly comprehensive review of the different approaches to examining 
these orbits. 

Hsiao and Scheeres [6] brought together concepts involving formation flying about a 
stable trajectory. Their work developed approaches to the implementation of distributed 
interferometric imaging by relying on the natural motion of spacecraft relative to each other. 
They were interested in the relative dynamics over a time period much less than the orbital 
period. Additionally, they defined preferred planes of motion for their stabilized system. 

While these analyses investigated spacecraft motion in the vicinity of a libration point, 
and even formation flying in that vicinity, they did not provide an analytical description 
of the relative motion of such spacecraft. Recently, NASA’s Goddard Space Flight center 
has begun developing several multi-spacecraft missions to the vicinity of the Sun- Earth L<i 
point. The planning for these missions requires just that type of relative motion description, 
suitable for orbit control design and analysis. The mission concept is that of a very sparse 
aperture plane located near the L 2 point, with a hub spacecraft located at the plane’s 
center. The hub is surrounded by free-flying telescope spacecraft in loosely a hub-and- 
spoke configuration. It is desired that the description of the telescope motion with respect 
to the hub be in a frame of reference which is amenable to the description of the mission’s 
scientific observation requirements. 

The present research is motivated by formation flying concepts for one particular such 
mission, the MicroArcsecond X-ray Imaging Mission (MAXIM) Pathfinder. MAXIM Path- 
finder consists of perhaps six telescope spacecraft forming a planar aperture, along with an 
optical hub located at the center of the configuration. The individual telescopes are located 
from 100 to 500 meters from the hub, and evenly spaced around the aperture plane. There 
is also a detector spacecraft located approximately 20,000 km and 90° out of the plane 
formed by the flat aperture. The entire system is to be in orbit about the Sun-Earth L 2 
point. 

In support of this mission, the circular restricted problem was used as a basis for 
the derivation of differential equations which describe the motion of a typical telescope 
spacecraft relative to the optical hub. Following the approach of Richardson, a modified 
Lindstedt-Poincare method was used to develop an uncontrolled periodic solution for this 
relative motion, which is a requirement for maintaining the telescope formation. The so- 
lution includes linear effects of the hub motion about L<i\ quadratic hub motion effects 
are used to develop relationships between the telescope frequency and amplitude. In the 
course of the analysis, a halo-type orbit of the telescopes about the hub is chosen to provide 
periodic motion in the aperture plane. Comparison is made between the resulting relative 
motion solution and a baseline numerical solution of the full circular restricted equations of 
motion. 

The selection of the initial state is one of the most difficult problems in the simulation. 
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As demonstrated by the analytical solution, arbitrary selection of the initial conditions 
generally leads to divergence of the trajectory. Even with careful selection of the initial 
state, the orbit might not close in the typical 200-day period; however, it will likely exhibit 
a useful trajectory over, say, 15 days. It can be assumed that as the aperture plane orbits 
the L 2 point, a typical astronomical observation will last only a few days, before the plane 
is reoriented to observe some other target. The application of orbit control, at this time, 
serves to reset initial conditions. 


DIFFERENTIAL EQUATIONS 

Consider a body in the vicinity of the Sun- Earth L 2 libration point, acted upon by the 
force model of the classical circular restricted three-body problem. In this model, depicted 
in Figure 1, the Earth-Moon barycenter is treated as being in a circular orbit about the 
Sun, the spacecraft mass is considered to be negligible as compared to the two primaries, 
and only point-mass gravitational forces are considered. 




For a system of two such spacecraft, a telescope and the optical hub, the differential 
equations of motion for the telescope relative to the hub are 


_ f T t r k \ ( T t T h \ 

\plt 3 plh 3 ) \P2t 3 P2h 3 ) 

- + D '> 4 ■ -■ miI ‘ - D1) (s? “ £?) *■ 


( 1 ) 


where 


r — vector from hub to telescope = r< — 
n — vector from L 2 to object i 
fix = solar Keplerian constant 

fi 2 = terrestrial Keplerian constant (Earth + Moon) 
pu — distance from Sun to object i 
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p 2 i = distance from Earth-Moon barycenter to object i 
x e = distance from system barycenter to L 2 
D\ = distance from system barycenter to Sun 
D 2 = distance from system barycenter to Earth-Moon barycenter 
x = unit vector parallel to Sun-Earth line of syzygy, (Sun-Earth direction). 

The subscripts h and t refer to the hub and telescope. 

The differential equations of motion are now expanded in powers of the distances between 
L 2 and the hub, as well as between the hub and the telescope. This is done with the intention 
of developing an analytical solution in terms of these quantities, useful for performing a 
control system analysis. Using binomial expansions, the telescope motion relative to the 
hub is given by 


r = -pi 


r h (p'e 

{x e + 


“ ^2 




(x e 




where 

A { rh \ 2 . 2 ^/» 

\x e + D\ ) x e + D\ 
. A r 2 + 2th • r 2x 
1 “ (x e + D\) 2 + x e + D! 


^ ( rh V 4. 

2 _ \xe-D 2 ) + x e -D 2 
A ^ — + 2r/ » ' r . 2j 

2 (i e - -D 2 ) 2 x e — D2 


( 2 ) 


A magnitude ordering system is now employed in order to truncate Eq. (2), prior to 
forming the solution. It is noted that the motion of the system takes place on two separate 
distance scales. There is a large distance scale, in which the motion of the hub about L 2 
is described, relative to the motion of L 2 within the context of the three-body problem. 
Then, there is a small distance scale, in which the motion of the telescope about the hub is 
described, relative to the motion of the hub about i 2 - 


Table 1 

PHYSICAL AND DERIVED CONSTANTS [7] 


VI 

132,712,440,017.987 km 

X e 

151,105,099.094445 km 

V2 

403,503.236 km 3 /s 2 

D 1 

454.84086785372 km 

n 

0.199106385xl0“ 6 rad/s 

D 2 

149,597,415.850132 km 
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Table 2 

DISTANCE RATIOS AND BASIC ACCELERATIONS 


r/r h 

8.333 x 10~ 7 

T h /{x e + Di) 

3.971 x NT 3 

r h /(x e - D 2 ) 

3.980 x 10" 1 

r/(x e + Di) 

3.309 x 10“ 9 

r/(x e - D 2 ) 

3.316 x 10~ 7 

^/(x e +D,) 2 

5.812 x 10" 6 km/s 2 

M/{x e -D 2 f 

1.775 x 10~ 7 km/s 2 


For ordering purposes, consider the distance between the hub and L<i to be on the 
order of 600,000 km, and the distance between the hub and telescope to be on the order of 
500 m. Using the constants of Table 1, where n is the terrestrial mean motion about the 
Sun (assumed constant), the relative distances are approximated by the ratios of Table 2. 
In the differential equations, these ratios scale the basic acceleration quantities which also 
appear in the table. Terms involving these basic accelerations scaled by rh/{x e + -Di) and 
Th/{^e — D 2 ) are now designated as being of order 1; terms involving the basic accelerations 
scaled by r/(x e 4 £>i) and r/{x e — D%) are designated as being of order 3. 

Retaining terms in Eq. (2) only through order 3, substantial algebra gives the truncated 
differential equations as 

t - A [— r 4- 3xx] 4 B [3xr* 4 3:^ 4 (3^ ■ r - ISxz/Jx] 

+ C[(3r h • r - 15 xx h )i h + § (r h 2 - 5x h 2 )r - ^{2x h T h • r - 7xx h 2 + ar>, 2 )x] , 

where 

A= a 1 1 & 

(x e + D i)3 + (s e -D 2 )3 

r, Pi | M2 

(Xe + Dtf + ( x e -D 2 )* 

r= m , a 2 

(^4Di) 5 (x e - JD 2 ) 5 ’ 

Note that this truncation includes terms which are linear in the coordinates of r and no 
more than quadratic in the coordinates of r*. Terms involving A, I?, and C are, respectively, 
of orders 3, 4, and 5; lower order terms do not appear. 

The acceleration vector r may be written relative to a rotating coordinate system which 
rotates at the constant angular rate n about the z-axis normal to the ecliptic, and with the 
x direction as previously defined. This gives 


where column vector notation is used to indicate the xyz vector components. 


x — 2 ny — n 2 x " 
y 4 2 nx — n 2 y , 
z 
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Consider now the linear motion of the telescope about the hub. From Eq. (3), the linear 
differential equations are given by 

it — A [— r + 3xx] , 


or, in component form, 


x — 2ny — (n 2 -f 2A)x 
y + 2 nx — (n 2 — A)y 
z + Az 


= 0 . 


( 4 ) 


As would be expected, these differential equations take exactly the same form as the 
linearized equations for the motion of a single object about L 2 - Clearly, the out-of-plane 
motion in the z direction is decoupled from the motion in the xy plane. The z equation 
describes simple harmonic motion, with a solution that may be written as 


z — A z sin (vt + VO? 


where v 2 = A. 

The fourth-degree characteristic equation for the in-plane motion is 

s 4 — (A — 2n 2 )s 2 — (n 2 4* 2 A) (A — n 2 ) — 0. (5) 

Keeping in mind the relative magnitudes of n (« 0.0172 rad/day) and A (« 0.0012 rad/day), 
the solution to this equation has four distinct roots. Two roots, positive and negative real, 
correspond respectively to a divergent and convergent mode. The remaining two roots are a 
purely imaginary conjugate pair, corresponding to oscillatory motion with natural frequency 
A given by 

A 2 — +n 2 + iJ(j — n 2 ) 2 + (n 2 + 2A)(A n 2 ). 

Periodic motion occurs if the in-plane initial conditions are selected so as to excite only 
the oscillatory linear modes. From the solution of the eigenvalue problem, placing this 
requirement upon the initial conditions gives 

x(0) — Xy(0)/k and y(0) = -k\x(0), 

where 

A 2 + n 2 + 2 A 
k ~ 2 An 

Under this relationship among the initial conditions, the complete linear solution for the 
telescope motion about the hub may now be written in the form 

x = — A x cos(A t + <f>) 
y = kA x sin(At + cf>) 
z — A z sin(^t -h -0)* 

It is desired that the constellation of telescope spacecraft remain in approximately the 
same specified plane over the few days’ duration required for observations. To achieve 
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this orientation, a halo-type orbit of the telescopes about the hub is selected. Such an orbit 
provides approximately periodic motion in the aperture plane, with the out-of-ecliptic-plane 
fundamental frequency equal to A, the fundamental frequency of the xy- planar motion. 

In order to do this, the inclusion of higher-order forces is used to adjust the out-of-plane 
fundamental frequency. Consider the 2-component of Eq. (3): 

z = -Az + B(3xz h + 3 zxh) 4- C(~l2xx h z h 4- 3 yyh*h ~ 6 zx h 2 + \zy 2 + \zz h 2 ). 

Recall that A = i/ 2 , and define A such that A = A 2 — v 2 . Then, this differential equation 
may be written as 

z 4* A 2 2 = A2 + B(3xz h + 3zxh) 4- C(-12xxhZ h 4- 3 yyhZh ~ 6 zx h 2 4- \zyh 4- \zz h 2 ). 

Here, the magnitude of A allows the term A2 to be treated as a higher order term, grouped 
with the terms containing the coefficient C. Using this formulation, the linear solution now 
becomes 


2 = A z sin(Xt 4- ^), 

with the same fundamental period as the in-plane motion. 

Together with the x and y components of Eq. (3), the system differential equations are 

x — 2 ny — n 2 x — —2Ax 4 - B(— 6xx^ 4 - 3 yyn 4 - 3zzh) 

+ C( Ylxxh ~ 12 yx h y h - 12 zx h z h - 6 xy h 2 - 6 xz h 2 ) 
y + 2nx — n 2 y = -Ay + B{3xy h + 3 yx h ) 

+ C(-12xx h y h + | yy h 2 + 3 zy h z h - 6 yx h 2 + \yz h 2 ) 
z + A 2 z = Az + B(3xzh + 3zxh) 

+ C(-I2xx h z k + Zyyhz h - 6 zx h 2 + \zy h 2 + | zz h 2 ). 

LINDSTEDT-POINCARE DEVELOPMENT 

The analytical solution to the expanded equations of motion of Eqs. (6) is now developed 
using a modified version of the Lindstedt-Poincare method. This method involves a se- 
quential solution of a system of differential equations, ordered by magnitude of the terms. 
Simultaneously, periodic motion is ensured through the placement of restrictions upon the 
initial conditions. In this development, the equations are expanded through third order in 
the small quantities, which is defined as terms that are at most linear in the motion of the 
telescope relative to the hub, and quadratic in the motion of the hub relative to £2. 

Introduce the non-linear frequency terms by series expansion, and change the indepen- 
dent variable from t to r, where r = u)t, and the asymptotic series 


u) — 1 4- u)\ 4- 4- * * * 


is used to scale the linear frequency. 

Recall that the forcing terms in the differential Eqs. (3) are ordered such that the 
terms involving A, B, and C are, respectively, of orders 3, 4, and 5. Now, for notational 
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convenience, these terms are reordered as orders 1,2, and 3. The solution vector is assumed 
to take the form of an asymptotic series, such that 


’ £ ' 


x\ 4 £2 4 £3 4 * • * 

v 

— 

V\ + 2/2 4- V3 4 — 

z 


Zi 4- Z2 4- z 3 4* * • • 


where the ordering by subscript is consistent with the reordering of the terms of the differ- 
ential equations; the order of a given term is specified by the subscript. 

This expansion is then substituted into the differential equations of Eqs. (6), and terms 
collected by order. Using primes to denote differentiation with respect to r, the resulting 
first-order equations for £i, yi, and z x , are given by 


" rr'i — 2 ny[ — (n 2 + 2A)x\ 
y" + 2 nx\ + (A - n 2 )y\ 
Z'{ + A 2 2 2 


= 0 . 


( 7 ) 


Note that these equations are identical in form to the linear equations of Eq. (4), with A 
replacing v in the z component equation. The second order equations contain contributions 
from the motion of the hub about I>2- As shown by Richardson, the hub motion relative to 
L 2 may also be expressed in a series form: 


Xh 


X\h + X 2 ft H 

Vh 

~ 

Vik 4- y2h 4 

. z h . 


z lh 4" Z2h 4 


In the second order equations, only x^, y\h, and z\h are included; at third order, the 2 h 
terms contribute as well. 

The second order equations for £2, y2? and 22, containing terms which are linear in the 
position of the hub, are then 

" x 2 - 2n^2 — (n 2 4* 2A)x2 
y*2 + 2n^2 + (A - n 2 )y 2 
z*2 + A 2 Z2 

( 8 ) 


—2u)i(x f { - ny[) + W{-2x x xi h 4- y x y xh 4- z x z xh ) 
-2wi (y x + nx[) 4- 4- y x x lh ) 

—2u)\z*{ 4- ZB{x x z xh 4- z x xih) 


Here, the telescope motion terms on the right side are now assumed to be known from the 
solution to the linear equations. 

Next, the third order equations for £3,2/3, and >2:3, containing terms linear in the 2 h hub 
position terms and quadratic in the 1 h terms, are 

£3 - 2 m /3 — (n 2 4- 2 ^ 4 )x 3 = — 2u)\(x*2 ~ ny 2 ) — 2u>2{x x — ny[) — u x x x 

+ 3B{y 2 yih + * 2 * 1 /, ~ 2x 2 x xh - 2x x x 2h + y x y 2h + z x z 2h ) 

+ §C[xi(1x\ h - y\ h - z \ h ) - 2y x x lh y Xh - 2z x x xh z xh ] 

(9a) 

y'l + 2 nx' 3 + (A- n 2 )y 3 = -2u)\(y" + nx' 2 ) - 2w 2 (y" + nx[) - w 2 y" 

+ W{x 2 y xh + y 2 x xb + x\y 2 h + y x x 2b ) 

+ 3 C [-4xix lh yih + yi {-2xl b + \y\ h + \z\ h ) + z x y xh z xh ] 

(9b) 
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Z% 4 = —2u)\Z2 “ (2cJ 2 4 4 A^i 

4- 3 B(x 2 z xh + z 2 x xh 4 Xiz 2 k 4- z x x 2h ) 

4 3(7 [-^ircuzu 4 yiyihZih 4 *i (-2^ 4 § 2 /?* 4 f^u)] • 

(9c) 

Again, the right side terms axe assumed to be known from the lower order solutions. 

Note that the left side terms at each order take an identical form. Therefore, a particular 
solution may formed for a general periodic forcing function. That solution may then be 
applied to each of the actual forcing terms, and the complete solution given by superposition. 
Higher order homogeneous solutions are not required, as they are identical in form to the 
terms of the linear solution, and thus are already present in the complete solution. 
Accordingly, consider the general case of the forced system, given by 

" x tf — 2 ny f — (n 2 4 2 A)x Acos(qr 4 0) 

y ff 4 2 nx* 4 {A — n 2 )y = B sin(yr 4 6) . ( 10 ) 

z" 4 \ 2 z _ _ T sin(pr 4 7 ) 

This vector is representative of all forcing terms which occur, where q and p are integers. 
Next, following the method of undetermined coefficients, the particular solution is assumed 
to take the form 

' x 1 I” R x cos {qr 4 6) 
y = Ry s\n(qr 4 9) . 

z J 1 Rz sin(pr 4 7) 

Substituting into the general system of Eq. (10), and solving for the solution amplitudes, 

_ 2 nqB — A{q 2 4 n 2 — A) 

Rx “ <z 4 4 {A - 2 n 2 )q 2 - (A - n 2 )(n 2 4 2 A) 

2nqA - B{q 2 4 n 2 4 2A) 

^ 4 (A - 2 n 2 )q 2 -(A- n 2 )(n 2 4 2 A) 

R * 
z ~ X 2 -p 2 ‘ 

Note that the denominator of R x and Ry is, in fact, the characteristic equation of the linear 
system, as seen in Eq. (5), where s 2 has been replaced by — q 2 . Therefore, upon examination 
of this particular solution, either q or p equaling A corresponds to the resonance condition. 
In that case, the forcing function has the same frequency as the homogeneous solution, and 
so the particular solution must instead have an amplitude which is secular in r. However, in 
order for the solution to be bounded for all r, such secular terms may not appear. Therefore, 
restrictions are placed on the solution, such that secular terms do not arise. 

For the case that p = A, the ^-component solution requires that T — 0. For q = A, the 
general x and y-component differential equations are combined into a single fourth-order 
differential equation: 

x ,tn — (A - 2 n 2 )x n — (A — n 2 )(n 2 4 2 A)x = [2 nBq 4 (A - n 2 - q 2 )A] cos (qr 4 9 ). 

It is then required that the forcing amplitude be zero, giving 

2 nBq 4 (A-n 2 - q 2 )A = 0 . ( 12 ) 
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Order 1 

The solution to the first order system of Eq. (7) is once again given by 


" Xi " 


" — A x cos(Ar + (f ) ) 

yi 

= 

kA x sin(Ar + </>) 

_ z i . 


A z sin(Ar 4- VO 


where here the solution is taken as a function of r. 


( 13 ) 


Order 2 


Here, in the second order equations, as previously mentioned, y\h, and z\h appear. 
Again from earlier analysis, the linear hub solution may be written as 




-A X h cos(A ht + 4>h) 

y\h 

— 

kA x h sin(A/,i + (j> h ) 

. z i/> . 


A zh sm(t/ h t + iph) 


Also, the hub frequencies A ^ and Vh may be written using a nonlinear asymptotic frequency 
scaling function u/j, where 


Uh — 1 + + UJ2h H • 

As already seen, the linear hub frequency is the same as that of the telescope; therefore 
A h ~ A Wh and For now, also assume that = u)\. 

Typically, a Lindstedt-Poincare analysis is employed for autonomous systems. Here, 
however, the forcing function is an explicit function of time, arising from the assumed 
solution for the hub motion. In the linear hub solution of Eq. (14), the independent variable 
t must now be written in terms of r, as 

t = t/lo 

- r/( 1 +u ) i +W 2 H ) 

== t (1 — 0^1 — 6^2 U ) I 2 + •••). 


Then, in the hub solution, 


A/j/ — 

= A(1 4* u>2h “* ^2 4- • ■ ■ )t. 

Define the second order contribution to this frequency correction as 


A 6^2 = U>2h ~ ^2, 


giving 


A ht — A(1 4- A(J2 + • ■ ■ )t. 

Similarly, for the out-of-plane frequency, 

Vht = Vh( 1 -W 1 -W 2 + wi 2 H )r 

= v h (\ + 
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where & represents the frequency correction through order 2. 

Also, write the hub phase angles (f>h and ^/i in terms of the telescope angles as 

(f) h = (f> -f A (j> and iph = 4- A 

where A <j> and A ip represent the offset in phase angles between the hub and telescope 
motion. Substituting into the linear hub solution of Eq. (14), and retaining the frequency 
correction through u )2 and u) 2 h, 


Xlh 


— A x h cos(A(l 4- Alo 2 )t + <f) + A$) 

yih 

— 

kA x h sin(A(l + Aaj 2 )t -I- (f) + A </>) 

- z l h . 


A Z h sin(^(l - u>)r 4* 4- A^) 


This linear hub solution and the first order telescope solution of Eq. (13) are substituted 
into the right side of the second order telescope equation, Eq. (8). After resolving the angles, 
the three component equations are then written as 

x% — 2ny 2 — (n 2 -I- 2A)x2 = — 2u>i A(A — nk)A x cos(Ar + <j>) 

— 3(l 4* 4>-) BA x A X h cos(A(2 4- Ao^t 4- 2(f) + A0) 

— 3(1 - y )BA x A x h cos(AAo;2T + A <p) 

+ ^BA z A Z h cos((A — Uh{\ -Cj))t - A ip) 

— %BA z A z hCos((\ + v>h(l —u)))t + 2ip + Aip) 

3/2 + 2 nx' 2 + {A — n 2 )j /2 = — 2uqA(— fcA + n)A x sin(Ar + <p) 

— 3BkA x A x h sin(A(2 + Aw? )r + 2<p + A <p) 
z" + A 2 Z 2 = 2ujiX 2 A z sin(Ar + ip) 

— \BA x A z h sin((A + j//,(1 — w))t + <p + ip + A ip) 

+ \BA x A z h sin((A — ^(1 - u>))r + <p-ip - A ip) 

— ^ BA z A x h sin(A(2 + Au >2 )t + ip + <p + A <p) 

+ \BA z A x h sin(AAa»2T - ip + <p + A <p). 

Within Eqs. (15), each term of the forcing function is treated in the manner of the 
general approach given in Eqs. (10) and (11). First, consider the potentially resonant terms 
in Eqs. (15a) and (15b). The condition of Eq. (12), to avoid resonance terms, gives 

— 2o>iA.A I [2n(— kX + n)A + (A — n 2 — A 2 )(A — nfc)] = 0. 

In general, this condition can only be satisfied for io\ — 0. Next, consider the resonant-type 
terms in Eq. (15c). Here, the condition to avoid resonance is 

2w\X?A z = 0 

and, again, this can only be satisfied for = 0. Accordingly, uq is now taken to be zero. 
Note that this requirement also gives u> = which will be used from this point forward. 


(15a) 

(15b) 

(15c) 
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The particular solution to Eqs. (15) is then built using the method of undetermined 
coefficients, as previously described. The resulting solution takes the form 

X 2 = A x A x h [p 2 i cos(A(2 + Aw 2 )t + 2<f> + A <j>) + P22 cos(AAa;2T + A<f>)] 

+ A z A zh [P23 cos ((A - V h (\ - w 2 ))r - Arp) 

+ p 2 iA z A zfl cos ((A + v h (\ - w 2 ))r + 2rp + A^)] 

V2 = A x A x h [<r 2 1 sin(A(2 + A w 2 )r + 2<p 4- A <p) + <x 22 sin(AAo; 2 T + A^>)] 

+ j4 z ^i,[<T 2 3sin((A - v h (l - w 2 ))r — Arp) + <r 24 sin((A + i/*(l - o> 2 ))r + 2 rp + AV>)] 
z 2 = A x A zh [k 2 i sin((A + i//,(l - w 2 ))r + cp + ip + Arp) 

+ K 22 sin((A — 1^(1 — w 2 ))r + <p — rp — A^>)] 

+ A z A x h [k 2 3 sin(A(2 + Aw 2 )t + rp + <p + A <p) + /c 24 sin(AAw 2 r - rp + <p + Acp)] , 

(16) 


in terms of p, cr, and k coefficient functions. 

Order 3 

The solution procedure at order 3 is much the same as for order 2, with much more compli- 
cated expressions. The previously determined lower order telescope motion is substituted 
into Eqs. (9), along with expressions for the hub motion. 

In the present analysis, the differential equations are examined only for potentially 
resonant terms — those terms in which A appears as the frequency in the r domain. Con- 
sideration of these terms will provide relationships between the acceptable values of A x and 
A z , as well as corresponding expressions for (Further research will allow the develop- 
ment of the actual order-3 solution terms.) The Cartesian components of these potentially 
resonant terms are denoted by 7Z X , lZ y , and 7^, where it is noted that Tl x and 7 Z y contain 
A x as a linear factor, while 7 Z z contains A z . 

Now examine the conditions to avoid resonance, given by Eq. (12), and the simpler 
z-component condition. It is clear that the coefficients of the resonant terms cannot satisfy 
these conditions with any degree of flexibility in the oscillation amplitudes — either or 
both of A x and A z would be required to be zero. Therefore, it is necessary to pursue an 
alternative. 

It is possible to increase the flexibility among the orbital parameters by introducing 
additional resonance-inducing terms. One way to do this is to require that the hub be in a 
halo orbit about 1,2- This means that — A/*, which gives 

Vh{ 1 ■“ ^2) — A(1 + Ao>2) 

(which may now also be used in the p, cr, and k functions). 

It can be seen that such an inclusion of additional terms does indeed introduce more 
flexibility. In particular, A x appears in the new z-component terms, while it is not in 1Z z ; 
similarly, A z is introduced into the x and y direction terms. 

Next, in order to allow these new terms to combine with the original resonance terms, 
it is required that the phase angles match. All in-plane trigonometric arguments are now 
required to be Ar + cf>; out-of-plane arguments are required to be At + i/>. From the hub halo 
motion analysis of [3], the hub phase angles must satisfy the relationship fa — fa = ^ r, 
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for arbitrary integer j. With this in mind, the trigonometric requirements lead to the 
relationship 

= (e + $)n, 

for arbitrary integer L 

Richardson also indicates that j must be odd (1 or 3) in order to avoid enormous hub 
motion amplitudes which would violate the series expansions being employed. Therefore, 
the same assumption is made here, that j must be odd. 

Using these results, 7^, 7 Z y , and 1Z Z are augmented with the additional resonant terms. 
For use in the equations at third order, the p, cr, and k functions need only be expanded 
through linear terms in Aa;2. Keeping this in mind, and recalling that Au;2 is defined as 
the difference between u)<ih and the augmented expressions may be written in the form 

1Z X — {<*laA x + + C*2a^4z 4 <^26^2^z) COS(Ar 4 <f>) 

TZ y ~ (PiaA x 4 Pib^ 2 ^x + /% aA z 4 fh.b^Az) sin(Ar + 4 >) 

H z = (71 aA x + 716^2^1 + 72a4z + 726^2^) (2 - j')(-l)* cos (At + <f>), 

where the a , /?, and 7 coefficients are functions of constants and the the hub solution 
parameters u) 2 a, A x h, and A z h- 

Recall that the conditions to avoid resonance are as specified by Eq. (12) and by the 
requirement that the z component of the resonant forcing function have zero amplitude. 
Here, these requirements now take the form 

2n(/3i a Ar 4 /3ibUJ2A x 4* $2 aA z 4 02b&2A z )\ 

4 (A — n 2 — )?){ot\ a A x + Otib(V 2 Ax 4 & 2 aA z 4* 0 ^ 2 b^ 2 A z ) = 0 (17) 

and 

l\aA x -I- 716^2^3; 4 72 a A z -I- 726^2^4z = 0. (18) 

Eq. (17) may be written in the form 

Cla^E 4- (lb^Ax + (2a A z 4 (2b^A z = 0. (19) 

Eqs. (18) and (19) provide two equations in the three variables A x , A Z1 and o?2- However, 
a third relationship may be introduced. Consider the order-1 solution of Eq. (13). Using 
the frequency and phase relationships previously discussed, 


X X m 


-A x cos(A r + <f>) 

yi 

— 

kA x sin(Ar + (f> ) 

*1 . 


. {2-j)(-l) e A z cos(\T + <f>) _ 


with j either 1 or 3, and £ either 0 or 1. This part of the solution corresponds to the case 
where the aperture plane maintains a roll angle ( about the y-axis. The range of this roll 
angle is related to the integers j and £ in the fashion shown in Table 3. 

Therefore, the ratio (2 — j)(—l) l A z /A x may be taken as an approximation to tan£. Let 
rj represent A z /A x . Assuming A x to be nonzero, Eqs. (18) and (19) may be written in terms 
of rj as 


7l a 4 716^2 4 72a ^7 + 726^2*7 = 0. 


( 20 ) 


13 



Table 3 

RELATIONSHIP OF APERTURE PLANE ROLL ANGLE TO j AND i 


3 

1 

* 

1 

0 

(0, 7t/2) 

1 

1 

(— tt/2, 0) 

3 

0 

(-• tt/2,0) 

3 

1 

(0, tt/2) 


and 


Cla + Cl6^2 4* £ 2 aV + C2b^2V = 0 (21) 

These two equations may be solved simultaneously for W 2 and 77. The combination of 
the equations leads to a quadratic, with two sets of solutions for L 02 and 77. Depending on 
the hub motion amplitudes, the solutions for rj correspond to approximate aperture plane 
roll angles from — vr/2 to 7 t/ 2, each associated with an infinite set of (A x , A z ) pairs. Each 
pair is, in turn, associated with a set of initial conditions for the relative telescope motion. 

A rule of thumb may be developed in order to bound the selection of A x and A z . 
Examining the solution to the linear telescope equations, given by Eq. (13), it can be seen 
that the approximate maximum distance of a telescope from the hub is 

dmax = A x max( y/l + 77 2 , k ) . 

Therefore, if the maximum distance is required to be no greater than a distance d max , A x 
may be selected such that 

A x < d m ax / max(v/l + t? 2 , k), 

and then A z = r]A x . 

Solution Summary 

The solution through order 2 is now constructed by adding Eqs. (13) and (16), and using 
the order-3 frequency and phase angle relationships. This combined solution is then 

x — —A x cos(Ar 4- 4>) 

+ ( P 2 \A x A X h — p 24 AzAzh(— 1)*) cos(A(2 + Aw 2 )r + <f>h + <f>) (22a) 

+ {pnA x A xh + P 23 A z A zh (-l) e ) cos(AAw2T + 4>h ~ 4>) 
y = kA x sin(Ar + <f>) 

+ {<72iA x A xh - (T 24 A z A zh (-l) e ) sin (A (2 + Alo 2 )t + <j>h + <£) (22b) 

+ {<722A x A x h - <J2zA z A z h{—l) 1 ) sin(AAw2T + (f>h~ <P) 
z = [(-1)^* cos(Ar + <f>) 

4- (^ 2 iA x A Z h + K23A z A x h(—lY) cos(A(2 + Au> 2 )r + <j>h + <f>) (22c) 

- (k 2 2 A x A z h + K 24 A z A xh (-l) 1 ) cos(AA u 2 t + <fih- <f>)\ (2 - j), 
where r = (1 + W 2 1). 
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SIMULATION 


The simulation results presented in this section compare the analytical relative motion 
solution with numerical integration. The motions of the hub and only one telescope were 
simulated to permit simpler interpretation. Given that the trajectory of a nonmaneuvering 
telescope will diverge from the nearby reference path about the hub and that finding initial 
conditions that produced useful trajectories was difficult, implementation of short duration 
flight times was fruitful. Assume that as the aperture plane orbits the L 2 point, a typical 
observation will last only a few days before reorienting the plane (and, if necessary, the 
orbit path) to observe some other target. By simulation of motion of 20 days duration, 
small diverging motions from the nominal halo orbit were maintained. 

The baseline solution is the numerical integration of the full equations of motion of the 
hub subtracted from that of the telescope, as shown in Eq. (1). In essence, there are two 
separate spacecraft orbiting the L 2 point with very similar initial conditions. The slight 
differences in initial conditions are due to the placement of the telescope with respect to the 
hub. The numerical calculations used to set up the other two solutions — the numerical 
solution to the truncated nonlinear equation of motion, Eqs. (3), and the analytical solution 
of Eqs. (22) — are outlined here. 

The algorithms are implemented as follows: 

• begin with a selection of A z ^ the amplitude of the motion along the z - axis of the hub 
with respect to L 2 

• use a consistent set of physical quantities taken from Dunham and Muhonen [7] for 
the constants used in computing A, £, C of Eqs. ( 3 ) and for A and k 

• compute the 7 and ( coefficients used in Eqs. (20) and (21) 

• solve for both of the 77 and u >2 solutions from Eqs. (20) and (21) 

• form two A x , A z = 77 A x pairs 


Table 4 

SUMMARY OF INPUTS AND RESULTING VALUES USED TO SELECT ORBIT 

SIZE AND ORIENTATION 


■A-xh 

(km) 

(km) 

3 

<jJ2 

U?2 + 

n 


r 

(deg) 

r 

(deg) 

A z 

(m) 

(km) 

280,667 

550,459 

1 

-0.29211 

-0.29185 

1.6115 

1.6147 

-58.18 

-58.23 

80.58 

0.08073 

280,667 

550,459 

3 

-0.29211 

-0.29185 

1.6115 

1.6147 

58.18 

58.23 

80.58 

0.08073 

227,219 

250,000 

1 

-0.26069 

-0.11043 

0.45009 

3.9851 

-24.23 

-75.92 

22.50 

0.1993 

227,219 

250,000 

3 

-0.26069 

-0.11043 

0.45009 

3.9851 

24.23 

75.92 

22.50 

0.1993 

211,126 

1,000 

1 

-0.23166 

-0.07373 

0.00178 

905.98 

-0.1020 

-89.94 

0.08902 

45.30 

211,126 

1,000 

3 

-0.23166 

-0.07373 

0.00178 

905.98 

0.1020 

89.94 

0.08902 

45.30 


Table 4 shows three combinations of amplitude selection, selections of j, and the results 
of calculations based on these selections. The application of Eq. (G-l) from Richardson [ 3 ] 
constrains the A X h corresponding to the selected A (This equation is plotted in Richard- 
son’s Figure H-l.) The amplitude of the hub motion along the 2 -axis was first chosen and 
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then the amplitude of the hub motion along the x-axis is constrained by this relationship, 
which assures halo orbit solution via correct frequency selection. The choice of A Z h of 
550,459 km is approximately the largest rvalue that can be selected and yield real roots 
for o> 2 * An intermediate value of 250,000 km and a very small value of 1,000 km were also 
selected for presentation. 

There are only two acceptable choices for j, which controls the sign on the z-position 
equation. The selection of j = 1 corresponds to a left roll orientation of the aperture plane; 
the selection of j = 3 corresponds to a right roll. This table shows the results of calculations 
based on fixing t = 1 and A x = 50m. (Results for l = 0 are not shown because the 77 values 
would be negative in Table 4 and lead to negative values of A z .) 

Because Eqs. (20) and (21) admit two solution pairs, these two pairs are each examined. 
The superscripts — and H- refer to the sign on the square root used in forming each solution. 
Two values are calculated for 77 and, in turn, the arctangent of 77, to yield the approximate 
aperture plane roll angle, 

The resulting magnitudes of the 77 terms, and their corresponding angle values should 
be highlighted. As suggested above, the large value of A z h is associated with the case in 
which the two roots for lo 2 approach a single double root. Accordingly, the resulting values 
of 77, £, and A z are also nearly identical. For the intermediate and small values of A z h , the 
quadratic gives distinct value pairs for u )2 and the approximate roll angle £. For the small 
Azh) the two resulting roll angles present extremes of the possible aperture plane orientation 
— based on 77”, the halo plane is nearly in the ecliptic ( x-y ) plane and based on 77*, the 
halo plane is nearly in the y-z plane. 

Because the user can select either A z h or the desired roll orientation, the user does not 
have control over the combination of halo plane size and orientation. When considering 
the roll orientation of the aperture plane, from — 7r/2 to tt/ 2, the user should realize that 
three months later, the orbit will have moved about the Sun such that the initial roll angle 
appears as a pitch rotation. This offers some orientation flexibility. Additionally, the user 
does not have complete freedom over the hub position and, of course, the aperture plane 
location with respect to £ 2 . 

Continue with the implementation of these algorithms: 

• compute the p, < 7 , k terms for Eq. (16) 

• compute the analytical expressions of Eqs. (22) 

The analytical solution is then formed by evaluating the analytical expressions for the 
motion of the telescope relative to the hub, as given by Eqs. (22). By differentiating these 
equations with respect to time, one can compute the full state. These equations may be 
used by trajectory control designers. 

For comparison purposes, the evaluation of this state at the initial time, t = 0, may be 
used to determine the initial state for use by two solutions using numerical integration: 

• Solution of Eq. (1) - full nonlinear equations of motion of the telescope with respect 
to the hub 

• Solution of Eq. (3) - truncated, to second order in the distance of the hub from £ 2 ? 
nonlinear equations of motion of the telescope with respect to the hub 
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Example Simulation and Results Comparison 

Presented here is an example that compares the results of the three solution methods for 
calculating the motion of the telescope with respect to the hub. 

“Pull” refers to the numerical integration of the separate equations of motion for the 
telescope and the hub, each spacecraft with respect to and then subtracting the hub 
motion from that of the telescope. This represents the implementation of Eq. (1). 

“Truncated” refers to the numerical integration of the truncated equations of motion 
for the telescope with respect to the hub. This represents the implementation of Eq. (3) 
with the terms that have a subscript h referring to the linear expressions for the motion of 
the hub. The linear hub solution is given by Eq. (14), which was further simplified to give a 
halo orbit setting both frequencies to the same value and setting both phase angles to zero. 

The initial conditions for both the full and truncated numerical integrations are obtained 
from the calculation, at t — 0, of the analytical solution. 

“Analytical” refers to the second order solution, which is an explicit function of time, 
given by Eqs. (22). 

Prom the intermediate case of Table 4, the out-of-plane linear hub amplitude from the 
1/2 point is selected to be 250,000 km ( A z h ). Using Eq. (G-l) of Richardson [3], the resulting 
value of A x k is 227,219.42 km. Then, the row with j — 1 or 3 is selected, corresponding 
to the desired aperture roll direction. The roll direction is implemented in the solution via 
Eq. (22c), and is also indicated in Table 3. Remembering that Table 4 was computed using 
l — 1, j = 1 is chosen for a negative roll or j = 3 for a positive roll. In this simulation, 
3 = 3- 

Select the initial amplitude between the hub and the telescope to be A x = 50 m, along 
the ar-axis. The A z ~ choice initially gives the telescope a linear amplitude of 22.5 m away 
from the hub along the z-axis. 

The analytical solution was first evaluated at t = 0 to obtain the initial conditions of 
the telescope relative to the hub, for use in the full and truncated solutions. This state is 
listed in Table 5, along with the initial conditions for the hub and phase angles. As needed 
for the numerical integration, the complete initial state for the telescope with respect to the 
1/2 point is then: 

x t (0) = £^(0) + z(0) ±t(0) = £/i(0) + i(0) 

!/t(0) = Vh{ 0) + 2/(0) yt{ 0) = yh( 0) 4- y(0) 

z t ( 0) = z h (Q) + *(0) z t (0) = z h (0) + z{ 0). 

Figures 2-4 show the small differences among the three solution types over a 20-day 
duration. The plots show the telescope motion with respect to the hub. Note that over 
5 days the solutions yield very similar results; the analytical solution is quite accurate for 
trajectory control analysis. It may be assumed that, as the aperture plane orbits the 
point, a typical observation will last only a few days, before reorienting the plane to observe 
some other target. By simulating runs of 20 days duration, one may see small diverging 
motions from the reference numerical solution. 

In Figure 2, observe that the full and truncated ^-component solutions are essentially 
the same. Considering the scale of the axis labeled “x position”, the analytical solution 
is close to the other solutions over 5 days. The curve labeled “Trunc C=0” shows the 
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Table 5 

INITIAL CONDITIONS OF TELESCOPE AND HUB 


*(0) 

-64.7796 m 

^(o) 

227,219.42 km 

2/(0) 

0.0 

m 

2//i(0) 

0.00 km 

*(0) 

26.4452 m 

2/1 (0) 

250,000.00 km 

i(0) 

0.0 

m/day 

ih (0) 

(X/k)y h (Q) km/day 

2/(0) 

4.4205 m/day 

2//i(0) 

-&Aa:/,(0) km/day 

m 

0.0 

m/day 

2 /. (0) 

-(2 - j)\A zh sm(<}> h ) km/day 

<f> - 

0.0 

rad 



4>h 

0.0 

rad 





Figure 2 Solution Comparison for x-Direction Telescope Motion 

numerical solution to the truncated equations, where the highest order terms are omitted. 
Clearly, for this example, the effects of these terms are not significant in the x-direction. 

In Figure 3, observe that the y-component solutions are all essentially the same for the 
scale shown. 

In Figure 4, it is seen that both the full and truncated ^-component solutions are quite 
close, with the analytical solution remaining close to these solutions over 5 days. Recall 
that the truncated solution has all of terms A, B, and C included. Here, note that the 
curve labeled “Trunc C=0” shows nearly the same result as the analytical solution. Now 
the importance of the higher order “C” terms is evident. It indicates that a higher order 
analytical solution could be expected to more accurately track the numerical solutions. 

SUMMARY AND CONCLUSION 

This report details preliminary work describing the formation flying between spacecraft 
near the Sun-Earth L 2 libration point, beginning with the circular restricted three-body 
problem for the hub motion about L 2 . 

The halo orbit was specified as a desirable aperture plane orbit, instead of a Lissajous 
orbit, because it is periodic. Additionally, over a one or two day observation period, the 
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Figure 3 Solution Comparison for y-Direction Telescope Motion 
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Figure 4 Solution Comparison for ^-Direction Telescope Motion 


spacecraft telescopes of the aperture plane will minimally separate due to natural motions. 

The analytical solution for the motion of a typical telescope relative to the hub is 
presented, including terms which are linear in the hub motion. These equations are for use 
in orbit control system design. This solution is linear in the position of the hub relative to 
L 2 ; quadratic hub position resonance effects are used to determine the relationships between 
the various parameters of the system orientation. 

Also developed was the full numerical solution of the relative telescope motion, as a 
reference trajectory. Additionally, the numerical solution to a truncated system of equa- 
tions for the telescope motion relative to the hub was formed for comparison. The three 
solutions were studied over a short 20 day duration, with the expectation that astronomical 
observations would be on the order of a few days, before reorienting the aperture plane. 

One example is presented and used to compare the accuracy of the analytical solution 
with the reference trajectory. Over 5 days, the analytical solution provides a reasonable 
approximation to that reference. This work provides a preliminary foothold for use in 
trajectory control design and analysis. This model is cast to show telescope motion relative 
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to a local hub-centered frame, rather than formation flying of multiple spacecraft relative 
to the libration point. 

Future work may involve the study of higher order contributions, such as terms which 
are quadratic in the hub motion about as well as deviations from the circular restricted 
problem. 
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